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Abstract 

The cores of edge dislocations, edge dislocation dipoles and edge dislocation loops in planar 
graphene have been studied by means of periodized discrete elasticity models. To build these 
models, we have found a way to discretize linear elasticity on a planar hexagonal lattice using com- 
binations of difference operators that do not involve symmetrically all the neighbors of an atom. 
At zero temperature, dynamically stable cores of edge dislocations may be heptagon-pentagon 
pairs (glide dislocations) or octagons (shuffle dislocations) depending on the choice of initial con- 
figuration. Possible cores of edge dislocation dipoles are vacancies, pentagon-octagon-pentagon 
divacancies, Stone- Wales defects and 7-5-5-7 defects. While symmetric vacancies, divacancies and 
7-5-5-7 defects are dynamically stable, asymmetric vacancies and 5-7-7-5 Stone- Wales defects seem 
to be unstable. 

PACS numbers: 61.72.Bb, 05.40.-a, 61.48.De 
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I. INTRODUCTION 



Graphenei^ and other two dimensional (2D) crystals' have been experimentally observed 
quite recently. This discovery has led to new physics where quantum relativistic phenomena 
can be mimicked and tested experimentally in condensed matter physics^. Among striking 
electronic properties of graphene due to its quantum electrodynamics-like spectrum, there 
are two chiral quantum Hall effects, minimum quantum conductivity in the limit of vanishing 
concentration of charge carriers and strong suppression of quantum interference effects. 
Ballistic transport on submicron distances at room temperature makes graphene a promising 
material for nanoelectronics 4 . 

Defects in graphene strongly affect its electronic and magnetic properties" 1 ^, which may 
be described by the Dirac equation on curved spaced Irradiation experiments show that 
pentagon-heptagon pairs (5-7 defects), vacancies, divacancies (5-8-5 defects comprising an 
octagon and two adjacent pentagons) and adatoms are commonly obtained, but seemingly 
not Stone- Wales (SW) defects (two adjacent 5-7 defects with opposite Burgers vectors and 
whose heptagons share one side, briefly 5-7-7-5 defects) 9 . The far field of 5-7 defects cor- 
responds to that of edge dislocations in elasticity, while the far field of vacancies and diva- 
cancies is that of an edge dislocation dipole. Quite recently, experimental observations of 
edge dislocations on high-quality graphene grown on Ir(lll) have been reported^. Studies 
of defects and their motion are important to assess the mechanical response of graphene at 
the atomic scale and, as indicated above, its electronic properties. 

One common way to describe defects in graphene is to use ab initio calculations. Density 
functional theory (DFT) has been used to ascertain the magnetic properties of graphene 
sheets and single-wall nanotubes with vacancies'. Local spin DFT has been used to describe 
glide and shuffle dislocations in irradiated graphitic structures^. Molecular dynamics (MD) 
has been used to discuss the stability of nanotubes under tension^ and also in the presence 
of different 5-7 pairs such as 5-7-7-5 and 7-5-5-7 (similar to SW but now the pentagons share 
one side) defects^. Atomistic Monte Carlo simulations are less costly than MD and have 
been used in studies of the stability of single graphene sheets^. Classical structural models 
of graphene may account for chirality effects in nanotubes and allow to assess the impact of 
the lattice structure on some elastic properties^. However, these classical models lack the 
ability to generate and move defects. The approach presented in this paper is different. 
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Relevant defects in graphene are the cores of different edge dislocation and edge dis- 
location dipoles. Thus we will regularize appropriately linear elasticity on the graphene 
honeycomb lattice and describe which are the stable cores of different edge dislocations and 
dipoles^. It is well known that cores of dislocations in crystals with covalent bonds are 
very narrow, so that the elastic field decays quite fast to that given by linear elasticity as 
the distance to the dislocation point (in 2D the dislocation line is a point) increases^ 1 ^. 
This means that we can regularize linear elasticity on a relatively small hexagonal lattice 
and impose boundary conditions corresponding to the elastic field of an edge dislocation (or 
a dislocation dipole) on a boundary which is sufficiently far from the dislocation point for 
the differences of the displacement vector to be well approximated by their corresponding 
differentials. The result will match seamlessly a calculation on a much larger lattice provided 
the far field of a dislocation is the same as that given by linear elasticity as we depart from 
the dislocation point. Despite the slow decay of the elastic strain away from the dislocation 
point, differences and differentials of the elastic displacement become indistinguishable a few 
atoms away from the dislocation point. 

Recently, we have developed periodized discrete elasticity models of dislocations in cubic 
crystals that describe their motion and interaction^^&^i. These models appear to provide 
the simplest correction to the equations of elasticity allowing nucleation and motion of defects 
and have two main ingredients. Firstly, by discretizing elasticity, a linear lattice model 
involving nearest and next-nearest neighbors is found. In the continuum limit, this lattice 
model provides the equations of elasticity with the appropiate crystal symmetry. Secondly, 
we need to account for the fact that certain atoms at the core of a dislocation change their 
next neighbors as the dislocation moves. Thus to describe dislocation motion we need an 
algorithm to update neighbors. Alternatively, we can periodize discrete differences along 
the primitive directions of the crystal by using an appropriate periodic function, thereby 
restoring crystal periodicity. This periodic function is selected in such a way that elasticity 
equations are recovered far from defect cores and stable static defects can be generated 
using their known elastic far fields at zero applied stress. For applied stresses surpassing the 
Peierls stress of the material, the defects should move and the value of the Peierls stress can 
be used to calibrate the periodic function- 
Here we extend this idea to the study of defects and their impact on the mechanical 
response of graphene sheets at zero temperature. Mostly planar graphene is very different 
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from cylindrical carbon nanotubes having very small radii and large curvature and from 3D 
graphite whose descriptions necessarily require studies different from the present one. Two 
facts complicate the task of designing lattice models of graphene. One is that graphene has 
a two-atom basis. The other is that a planar hexagonal lattice is intrinsically anisotropic 
even if its continuum limit is isotropic elasticity. We have dealt with these issues by using 
difference operators whose continuum limits are linearly independent combinations of the 
partial derivatives entering the Navier equations of linear elasticity. The key idea to choose 
our difference operators is that they do not have to involve all neighbors of an atom on an 
equal footing. Once we have an appropriate discretization of linear elasticity, we periodize 
the resulting lattice model in a way that allows dislocation gliding. Adding thermal effects 
and local curvature effects due to ripples^ 1 ^ increases the complexity of models and we have 
omitted these effects in the present paper, which is organized as follows. Section [III recalls 
some basic details on the structure of graphene lattices and presents the stable defects we 
have obtained by solving the periodized discrete elasticity models. The basis of these models 
are lattice models obtained by discretizing elasticity, as indicated in Section IIII1 Periodized 
discrete elasticity is discussed in Section HVJ Section [V] describes numerical tests of defects 
in graphene carried out with the full 2D model and with a scalar reduced version thereof. 
Among solutions of our models, we have found a stable octagon defect. Known defects 
such as pentagon-heptagon pairs and 5-7-7-5 SW defects move and interact as expected. In 
particular, the two 5-7 defects comprising the two edge dislocations with opposite Burgers 
vectors of a SW glide to each other on their common gliding line and annihilate. This is 
not the case if we have a 7-5-5-7 defect (similar to the SW, but now the pentagons share 
a common side) 1 -, which is stable because the dislocation centers of the component edge 
dislocations are displaced one atomic distance from each other in different glide lines^ 3 . As 
observed in experiments, we have also obtained stable 5-8-5 divacancies^ and symmetric 
vacancies^. We find that a symmetry-breaking vacancy^ evolves toward a simple threefold 
symmetric vacancy^. Lastly, Section [VI] contains our conclusions. 

II. STRUCTURAL CHARACTERISTICS OF GRAPHENE AND DEFECTS 

Figure [Tj illustrates the structure of a graphene sheet comprising a single-layer of graphite. 
This 2D hexagonal lattice is equivalent to a cubic lattice with a two-atom basis generated 
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FIG. 1: (Color online) Graphene hexagonal lattice. The unit cell vectors are a = a(l,0) and 



carbon atoms A = (0,0) and B = a(cos(^), sin(^)) belonging to two sublattices. An atom A has 
three nearest neighbors, six next-nearest neighbors and three second-nearest neighbors. 



FIG. 2: (Color online) Numerically generated single pentagon-heptagon defect and its motion. The 
pentagon marks the lower end of an extra column of hexagons and its basis forms the top of the 
heptagon, (a) Under sufficiently large applied shear stress, a dark atom moves to the right and 
a light one to the left as indicated by the arrows, (b) The result is that the pentagon- heptagon 
defect moves one step to the right. 

by two non orthogonal unit cell vectors, a = (l,0)a and b = (l,^/3)a/2, where a is the 
lattice constant. The length of a hexagon side is / = a/\^3. Dark and light colors are 
used to distinguish the two sublattices: dark atoms belong to sublattice 1 and light ones to 
sublattice 2. 

Graphene layers often contain defects: an experimental study of defects generated by 



b = a(cos(^), sin(^)), with a lattice constant a = 2.461 angstroms. The unit cell contains two 
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FIG. 3: (Color online) Numerically generated shuffle dislocation with dangling bond. An octagon 
marks the position of the core, (a) Under sufficiently large applied shear stress, a light atom moves 
to the right and a dark one to the left as indicated by the arrows, (b) The result is that the edge 
dislocation moves one step to the right. 

r~| 

irradiation and images thereof is reported in Ref. |9j. Several defects have been obtained by 
using the periodized discrete vectorial model presented in the next Section. They appear in 
Figures [21 El IH and In these figures and successive ones in the paper, only the positions of 
the atoms have been calculated numerically. The bonds between neighboring atoms are only 
aids to visualize the structure of the lattice. 

A simple defect in a graphene sheet is the single pentagon-heptagon (5-7) pair, depicted 
in Figures EJ^a) and (b). As we will show later, this defect represents an edge dislocation 
and its gliding motion can be characterized by atom motion together with breakup and 
attachment of bonds between atoms. Note that link breakup and union occurs only between 
atoms in the direction of defect motion. As a result of this motion, the neighbors of the 
moving atoms in the rows immediately above and below them change. If we recall that 
bonds between neighboring atoms are only visualization aids, saying that "one defect moves 
by breaking bonds between some atoms and creating bonds with others" is only a figure of 
speech, not a consequence of the model. 

Our models show that edge dislocations can have cores different from those depicted 
in Figure [2j In Fig. [31 eight atoms, one with a dangling bond, form the core of an edge 
dislocation. When the dark atom moves to the left as indicated by the arrow and is attached 
to the light atom with the dangling bond, the light atom formerly connected to it moves to 
the right and one of its bonds is left dangling. The overall result is the one-step motion of 
the octagon defect to the right. 

Figure H] depicts a vacancy defect which is a possible core of an edge dislocation dipole, 
as it will be shown later. Fig. H^a) is the initial configuration obtained by discretization of 
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FIG. 4: (Color online) Numerically generated edge dislocation dipoles. Arrangement of atoms in 
a vacancy: (a) initial configuration as given by linear elasticity; (b) final configuration after time 
relaxation in the unstressed lattice. 



FIG. 5: (Color online) Numerically generated edge dislocation dipoles. (a) 5-8-5 divacancy. (b) 
5-7-7-5 Stone- Wales defect. 

the elastic field of an edge dislocation dipole. Fig. H](b) is the final configuration obtained 
by evolution of (a). 

Figure 0(a) depicts another possible core of an edge dislocation dipole: a divacancy 
formed by an octagon with two adjacent pentagons. This 5-8-5 defect is dynamically stable. 
Yet another different core of the dipole is the Stone- Wales (SW) defect formed by a pair of 
heptagon-pentagons, as depicted in Fig. [5](b). This 5-7-7-5 defect is created by the SW bond 
rotation of one light atom forming the lowest vertex of an hexagon and bond reattachment 
to another light atom to form the top of a pentagon and, at the same time, the basis 
of an heptagon. This bond rotation leaves an oppositely oriented 5-7 defect next to the 
other one. This SW defect has been obtained by superposition of the initial guesses for the 
elastic displacement fields of two opposite edge dislocations with heptagon-pentagon cores 
that share the same glide line. In elasticity, we would expect that the two component edge 
dislocations of this dipole glide towards each other and annihilate, leaving the undisturbed 
lattice as a result. This is in fact what the numerical solution of our model shows. Thus the 
SW defect represents an edge dislocation dipole and it seems to be dynamically unstable^. 
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FIG. 6: (Color online) The two component 5-7 edge dislocations of a SW defect separate under a 
large applied shear stress. This figure has been calculated using a simple periodized scalar model 
in a lattice with 18 x 18 lattice spacings. 



If a sufficiently large shear stress is applied to a lattice that contains one SW defect, its two 
component dislocations drift apart as shown in Fig. O 

If we form a different dipole with component dislocations having opposite Burgers vectors 
and different glide lines, the two component edge dislocations glide towards each other but 
do not meet. Instead they form a stable dislocation dipole. As in the case of nanotubes^, 
the simplest realization of this idea is a 7-5-5-7 defect in which the two pentagons (not the 
two heptagons as in the SW configuration) share a common side. Our results show that the 
7-5-5-7 defect is dynamically stable^. 



III. DISCRETE ELASTICITY MODELS FOR GRAPHENE SHEETS 

We want to produce a lattice model that reduces to the equations of linear isotropic 
elasticity in the far field of a defect. In the continuum limit, elastic deformations of graphene 
sheets are described by the Navier equations for the 2D displacement vector (u,v): 

d 2 u _ d 2 u d 2 u d 2 v 
d 2 v _ d 2 v d 2 v d 2 u 

P W ~ Cm w + Cu W 2 + ( 66 + 12) cW (2) 

where p is the mass density. In the basal plane, graphite is isotropic, so that C^e = p, 
C\2 = A and C\\ = A + 2p, in which A and p are the Lame coefficients. 
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A. Discretization of the Navier equations and linear lattice model 



We would like to find a linear lattice model by discretizing the Navier equations (j2J). 
This is not an obvious task because any symmetric combination of differences involving 
either all nearest neighbors, or all nearest and next nearest neighbors, etc. only yields a 
multiple of the Laplacian in the continuum limit, no matter how many neighbors we use. We 
have allowed combinations of differences that do not involve all neighbors of the same type 
symmetrically. The idea is to select three difference operators that yield three independent 
linear combinations of d 2 u/dx 2 , d 2 u/dy 2 and d 2 u/dxdy in the continuum limit. Then we 
replace the partial derivatives of u in ([2]) by the combinations of our differences that provide 
those partial derivatives in the continuum limit. 

Let us select assign the coordinates (x, y) to the atom A in sublattice 1 (see Figured]). 
The three nearest neighbors of A belong to sublattice 2 and their cartesian coordinates are 
rii, n 2 and n% below. Its six next-nearest neighbors belong to sublattice 1 and their cartesian 
coordinates are rii, i — 4, . . . , 9: 




(3) 

Four of these atoms are separated from A by the primitive vectors ±a and ±b (see Figure 

H). 

Let us define four operators acting on functions of the coordinates (x, y) of node A: 

Tu = [u{n x ) - u{A)) + \u(n 2 ) - u{A)\ + [u{n 3 ) - u{A)), (4) 

Hu = [u{n 6 ) - u{A)) + \u{n 7 ) - u{A)), (5) 

D x u = [u(n 4 ) - u(A)\ + [u(n g ) - u(A)\, (6) 

D 2 u = [u{n 5 ) - u{A)\ + [u{n 8 ) - u{A)}. (7) 

Note that the operator T involves finite differences with the three next neighbors of A which 
belong to sublattice 2, whereas H and D\ involve differences between atoms belonging to 
the same sublattice along the primitive directions a and b, respectively. See Figure [3 D 2 



FIG. 7: (Color online) Neighbors of a given atom A. Only the neighbors labeled 1, 2, 3, 4, 6, 7 
and 9 are affected by the difference operators T, H and D% used in our discrete elasticity model 
with two slip directions. 



involves differences between atoms belonging to the same sublattice along a different choice 
of the basis vectors: a and c (parallel to the line joining atoms and n$ in Fig. 0, which 
is also another primitive direction). The operator D2 will be important when we want to 
consider a dislocation motion along slip directions a and c. Taylor expansions of these finite 
difference combinations about (x, y) yield 



Tu 




Hu 



d 2 u 




dx 2 




as a 



0. 



1. Model with three slip directions 



Let us assume that we want to allow dislocations to slip along any of the three primitive 
directions a, b or c. Then we replace in Equations Q and (T5J) Hu/a 2 , (4T — H)u/a 2 and 
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(Di — D 2 )u/ (\^3a 2 ) instead of d 2 u/dx 2 , d 2 u/dy 2 and d 2 u/dxdy, respectively, with similar 
substitutions for the derivatives of v. In terms of the Lame coefficients, we obtain the 
following equations at each point of the lattice: 

pa 2 ^ = 4// Tu + (A + fi) Hu + ^±=* (D 1 - D 2 )v, (8) 
pa 2 ^ = ApTv + (A + p) (AT - H)v + (D 1 - D 2 )u. (9) 

2. Model with two slip directions 

If we only allow slip along the two primitive directions a and b that form our vector 
basis, we should replace in Equations and (T5]) Hu/a 2 , (AT — H)u/a 2 and 2(Di — 3T + 
H i '2)u I '(V3a 2 ) instead of d 2 u/dx 2 , d 2 u/dy 2 and d 2 u/dxdy, respectively, with similar sub- 
stitutions for the derivatives of v. In terms of the Lame coefficients, the following equations 
are then obtained 

pa 2 ^ = ApTu + (A + p) Hu + U x - 3T + \h\ v, (10) 



pa 2 — = Ap Tv + (A + p) [AT - H)v + VA ^_ 7 (A - 3T + -#J u, (11) 

at every point of the lattice. To have slip directions a and c, we replace the operator 
-(D 2 -3T + H/2) instead of (D x - 3T + H/2) in dTOj) and (TU]). The same equations are 
found if B = (x,y) is an atom in sublattice 2. In this case, the relevant neighbors of B 
entering the definitions of T, H, D\ and D 2 have coordinates 

/ a 

n i = ( x ~ tt> 2/ + ) > n 2 = he + -, y + 





r» 2V3 

a 

it i = | '• - -■ I) — ! • ''<■<■ \ x +2 ,y ~ 

n 7 = (x + a,y), n 8 = [x-^,y+ ^y^j , n 9 = [ x + ~, y + ) . (12) 

We may see the hexagonal lattice as a set of atoms connected by springs. These springs 
connect each atom A with its nearest neighbors n%, n 2 and 713, and with its nearest neighbors 
along the primitive directions, n^, rie, m and rig. If we add symmetrically the missing two 
next-nearest neighbors and n 8 in the operator D\, its Taylor expansion produces the 2D 
Laplacian. Similarly, adding all the second-nearest neighbors, the Laplacian is found again: 
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Symmetric choices of neighbors only generate Laplacians, but not the terms involving cross 
derivatives. It seems reasonable to break the central symmetry about a given atom when 
defining finite differences by giving preference to the primitive directions. Notice that if we 
move along the lattice, the hexagonal arrangement itself is a source of anisotropy. Along the 
x direction (a primitive direction of the lattice), atoms are arranged in a 'zig-zag' pattern. 
The same arrangement occurs along the other primitive directions. However, along the y 
direction atoms are arranged in an 'arm-chair' pattern. 



B. Lattice model in primitive coordinates 

Equations fTTOT) and (TTTT) can be written in primitive coordinates u\, i = 1, 2 (with u[ = u', 
u 2 = v ')i by means of the transformation Ui = %jUj (summation over repeated indexes is 
intended), with 



(13) 



Writing Equations (JED and (J9j) or ( flOl) and ( TTTT) as pa d Ui/dt = LijUj, the linear equations 
of motion in the primitive coordinates are pa 2 d 2 u' i /dt 2 = L'^u'j with L'^ = T^T 1 L kn T n j, or 
equivalently: 

pa 2 ^ = ApTu' + (3H - D 1 + D 2 )u + (A + p) (h + Dl ~ ° 2 - 2t) v', (14) 





dt 2 ^ 3 v ' K n ' \ 3 

P" 2 ^ = ^ (A - D * ~ 3ffV + 4(A + 2p)Tv' + 2 -{\ + fi{D x - D 2 )u', (15) 

when there are three slip directions, or 
pa 2 d 2 u' A + p, 



2 dt 2 3 
pa 2 d 2 v' A + p 

2 dt 2 = 3 



[(if - Di)u' + (2/f + D a )i/] + T[(A + 3^)u' - 2(A + (16) 
[(H + 2Di)u' + (Di - H)v'\ + T[(A + 3p)v' - 2(A + //)«'], (17) 



when there are only two slip directions along the basis vectors a and b. Note that u' = 
(u — v/\f?>)/a and v' = 2v/(ay3) are nondimensional. Equations (fl4|) and (fT5|) do not look 
symmetric in the same way as ( fTBl and (fTTfl do because we have selected a basis along 
primitive directions a and b, so that c (which defines the third slip direction associated to 
the operator D 2 ) is not a basis vector. 
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C. Scalar model 



When the displacements in the y direction are negligible we may ignore the vertical 
component and work in cartesian coordinates. The evolution equations for the displacements 
in the horizontal direction are: 

d 2 u 

pa 2 — = ApTu+{\ + p)Hu. (18) 

This equation can also be obtained from u — (u 1 + v' /2) a by adding twice <HM to (fT7j) . and 
then setting v' = 0. The only slip direction of the scalar model is along a. 



IV. PERIODIZED DISCRETE ELASTICITY MODELS FOR GRAPHENE 
SHEETS 

The models described by (jHJ) - ©, ffTUl) - ffTTl) or (fl~8]) are linear and do not allow for 
the changes of neighbors involved in defect motion. An obvious way to achieve this is to 
update neighbors as a defect moves. Models such as (jHj) - (jSj), ( JTOT) - ( TTTT) or ( flBl) would 
have the same appearance, but the neighbors rii would be given by ([3]) and (TT2"]) only at the 
start. At each time step, we keep track of the position of the different atoms and update the 
coordinates of the rij. This is commonly done in Molecular Dynamics, as computations are 
actually carried out with only a certain number of neighbors. Convenient as updating is, its 
computational cost is high and analytical studies thereof are not easy. Another advantage 
of periodized discrete elasticity is that boundary conditions can be controlled efficiently to 
avoid spurious numerical reflections at boundaries. 

In simple geometries, we can avoid updating by introducing a periodic function of dif- 
ferences in the primitive directions that automatically describes link breakup and union 
associated with defect motion. Besides greatly reducing computational cost, the result- 
ing periodized discrete elasticity models allow analytical studies of defect depinning^ 1 ^ 1 ^, 
motion and nucleation^I. 

To restore crystal periodicity, we replace the linear operators T, H, Di and D 2 in (tHJ) - 
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(|T5|) or (|T6|) - (|T7I) by their periodic versions: 

T p u' = g{u'{n x ) - u\A)) + g(u'(n 2 ) - u'(A)) + 0( u '(n 3 ) - u'(A)), (19) 

= g(u'(n 6 ) - + ^(u'(n 7 ) - u'(A)), (20) 

D lp «' = g(u(n 4 ) - u\A)) + <?(«'(n 9 ) - u'(A)), (21) 

D 2pM ' = g(u(n 5 ) - u'(A)) + g(u'(n 8 ) - u\A)), (22) 

where (7 is a periodic function, with period one, and such that g(x) ~ x as x — > 0. In our 
tests we have taken g to be a periodic piecewise linear continuous function: 

. x, —a < x < a, 

g a (x) = { ' " " (23) 

2o -x + -V, a < x < 1 - a. 



l-2a 1 1-2q' 

The parameter a controls defect stability and mobility under applied stress, a should be 
sufficiently large for elementary defects (dislocations, vacancies) to be stable at zero applied 
stress, and sufficiently small for dislocations to move under reasonable applied stress^. The 
periodic function g can be replaced by a different type of periodic function to achieve a 
better fit to available experimental or numerical data. 

Periodized discrete elasticity is a Lagrangian model: the atoms are labeled from the start 
and we track their motion. The periodic functions allow us to simulate dislocation motion 
without updating neighbors thereby greatly reducing computational cost. 

In the simpler case of scalar elasticity ffl8|) . the corresponding periodized discrete elasticity 
model is: 

pa 2 — = 4pT p u + ( A + p)Hu, (24) 

( u{ ni )-u{A) \ ( u{n 2 )-u{A) \ ( u{n z )-u{A) \ 
T p u = ag^ \ + ag^ )+ag^ j. (25) 

The nonlinear function g is only needed for differences between the neighbors that may 
change due to defect motion. In the direction x, neighbors never change, therefore we use 
the operator H of Eq. (jSJ). Horizontal rows may shift, resulting in a shift of neighbors that 
is taken into account by using T p as in (125]) with the periodic function g a defined in (I23~l) . 



V. DEFECTS IN GRAPHENE 



In this Section, we discuss the defects obtained with our periodized discrete elasticity 
models. In our numerical calculations, we use for graphene the elastic constants of graphite 
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in the basal plain (which is isotropic), Cxi = C12 + 2C*66 = 1060 GPa, C12 = A = 180 GPa, 
Cm = fi = 440 GPa. 26 This yields v = 0.17. We have used a = 0.2 in the periodic function 
g a (x). For this value, the Peierls stress for a 5-7 defect is 0.025 /1, which is of the order 
of known values for covalent crystals^. Larger (smaller) values of a yield larger (smaller) 
Peierls stresses: for a between 0.15 and 0.3, the Peierls stress varies in the range between 
10~ 3 /^ and 10 _1 /i. We have not calibrated our model to a precise value but this can be done 
when Peierls stresses for single graphene sheets are measured in experiments. All defects 
have been calculated using the vectorial model (ITBT) - (fTTj) . periodized by means of (jl~9l) - (I2TI) . 
Analogous results can be found using appropriate computationally cheaper scalar models 
such as (1241) - (1251) . Scalar models are more convenient to analyze defect motion, interaction 
and even nucleation of defects constrained to move along a given primitive direction. In 
more complex geometries a combination of periodized discrete elasticity models and neighbor 
updating could be useful. 

In all cases, the construction of defects is similar. We use the displacement field of a given 
defect in continuum linear elasticity at zero applied stress both as initial and boundary 
condition for the discrete model. Then the model with overdamped dynamics (replacing 
second order time derivatives in the model equations by first order ones) is used to relax the 
initial field to the stable stationary solutions representing the sought defects. Applied stress 
can be implemented via the boundary conditions. To study defect motion and interaction, 
we can use the models with inertia. 

How large should our computational lattice be? We can have an idea by using results by 
Zhang et al2£. As the far field of defects, they considered linear elasticity with a strain energy 
that was a quadratic functional of the strain tensor and also of its first and second derivatives. 
They compared results given by this theory with those of classical linear elasticity and with 
results from atomistic simulations. In their figure 11, it can be seen that four lattice spacings 
away from the core of a SW defect or of a di vacancy in graphene, linear elasticity already 
approximates very well MD results (they use only 132 carbon atoms in their simulations). 
Thus a lattice with 18 x 18 lattice spacings (and therefore with 36 x 36 = 1296 atoms) 
should be sufficiently large to obtain good results for a centrally located defect by using 
our discrete elasticity models. The figures presented in this paper have been obtained using 
such a lattice and we have checked that our results do not change significantly by decreasing 
slightly the lattice size or by increasing it. 
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As we have said before, our model is more efficient than MD. In MD, one has two options 
when computing the force acting on each atom, either to include the interactions of one 
atom with all its neighbors, or to include only the neighbors within a certain distance of 
the atom and update them as time evolves. The cost of computing the force on each atom 
at each time step according to the first option is huge due to the large number of atoms 
involved, and therefore we need a very long time to compute just a few time steps. The 
second option (few neighbors and updating) is less computationally costly than the first one, 
but its cost is still large. In contrast to this, our model involves a few nearest neighbors and 
the periodic function used in the governing equations avoids the need of updating neighbors 
at each step. Thus computation of the forces at each step is cheap. We may compute the 
evolution of atoms in much larger lattices in a much shorter time. Since our forces are 
cheaper to compute, we may use higher order solvers to integrate the resulting system of 
differential equations. Since we are able to handle larger lattices, our trouble with numerical 
artifacts due to reflections on the walls of the lattice is smaller. In MD, reflections of waves 
at the boundaries of the computational domain limit the time spans over which simulations 
are reliable and there is a trade off between lattice size and the total time of a simulation. 
In fact, it is not yet clear how to prevent spurious reflections in MD simulations in an 
efficient manner-^. The simplicity of our model equations allows us to introduce simple non 
reflecting boundary conditions, thereby suppressing numerical size effects due to reflections 
at the boundaries. 

A. Edge dislocations 

Static edge dislocations can be generated using the overdamped version of (|T6l) -(fTTl) peri- 
odized by means of ( |T9i) - ( J2T1) and the elastic field of edge dislocations for (flOl)-(|TTT). To find 
the stationary edge dislocation at zero stress, we first write the corresponding stationary edge 
dislocation of isotropic continuum elasticity. The displacement vector u = (u(x, y),v(x, y), 0) 
of an edge dislocation directed along the z axis (dislocation line), and having Burgers vector 



(6, 0, 0) is 




(26) 
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cf. Ref. 
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pag. 57. (|26|) has a singularity oc r _1 at the core and it satisfies J c ((ix • V)u = 
— (5, 0, 0), for any closed curve C encircling the z axis. It is a solution of the planar stationary 
Navier equations with a singular source term: 

Au + y^7;V(V ■ u) = -(0, b, 0) 6(r). (27) 
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Here r = ^ ' x 2 + y 2 and v = A/ [2 (A + fi)] is the Poisson ratio; cf. page 114 of Ref. 

The continuum displacement ( 1261) yields the nondimensional static displacement vector 
in primitive coordinates U'(l,m) = [u(x,y) — v (x, y)/V3]/a, V'(l,m) = 2v(x : y)/(a\^3), 
where x = a(x' + y' /2), y = a^/3y' /2. The primitive coordinates x' = x' + I, y' = y' + m 
are centered in an appropriate point (x ; ,y' ) which is different from the origin to avoid the 
singularity in (!26j) to coincide with a lattice point. 

U'(Z,m) will be used to find the stationary edge dislocation of the discrete equations of 
motion. To this end, we replace the inertial terms pa 2 d 2 u' j 'dt 2 and pa 2 d 2 v' /dt 2 in (jlfip - 
( JT71) by pdu'/dt and fidv'/dt, respectively. The resulting overdamped equations have the 
same stationary solutions. We use an initial condition u'(7,m;0) = XJ'(l,m) given by Eqs. 
( j26|) . and boundary conditions u'(l,m;t) = U'(/,m) + (Fm,0) at the lattice boundaries 
(F is a dimensionless applied shear stress). If \F\ < F cs (F cs is the static Peierls stress for 
edge dislocations), the solution relaxes to a static edge dislocation (u'(l,m),v'(l,m)) with 
the appropriate continuum far field. 

Depending on the location of the singularity (x' , y' Q ), there are two possible configurations 
corresponding to the same edge dislocation in the continuum limit. Figure [2] shows the 
structure adopted by the deformed lattice (I + u'(l,m),m + v'(l,m)) when the singularity 
is placed between two atoms that form any non-vertical side of a given hexagon. The core 
of the dislocation is a 5-7 defect. If the singularity is placed in any other location different 
from a lattice point, the core of the singularity forms an octagon having one atom with 
a dangling bond, as shown in Fig. [3j The dangling bond causes this configuration to be 
more reactive due to the possibility of attaching impurity atoms to the dangling bond. The 
octagon can also be seen as a 5-7 defect with an extra atom inserted between heptagon and 
pentagon. The basin of attraction of the octagon configuration seems to be larger than that 
of the 5-7 defect. That the same dislocation type may have two different cores is a familiar 
fact in crystals with diamond structure and covalent bonds, such as silicon; see page 376 in 
Ref. ilTJ. There it is shown that the 60° edge dislocations may belong to the 'glide set' or to 
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the 'shuffle set'. Seen from a certain direction, the cores of the glide set look like 5-7 defects 
whereas the cores of the shuffle set look like octagons with a dangling bond attached to one 
of their atoms. Ewels et al use the names "glide dislocations" and "shuffle dislocations" for 
the 5-7 defects and the octagons with a dangling bond (glide + adatom), respectively^. 

The glide motion of edge dislocations occurs in the direction of their Burgers vector and on 
the glide plane defined by the Burgers vector and the dislocation line. In the configurations 
of Figs. [2] and [3], a supercritical applied shear stress will move the dislocation in the x 
direction on the glide plane xz. Our simulations show that the shuffle dislocations move 
more easily than the glide dislocations, as predicted by Ewels et al^. For conservative or 
damped dynamics, the applied shear stress has to surpass the static Peierls stress to depin a 
static dislocation, and a moving dislocation propagates provided the applied stress is larger 
than the dynamic Peierls stress (smaller than the static one) 22 . A moving dislocation is 
a discrete traveling wave advancing along the x axis, and having far field (u'(l — ct,m) + 
Fm, v'(l — ct, m.)). The analysis of depinning and motion of planar edge dislocations follows 



that explained in Ref. |22| with technical complications due to our more complex discrete 
model. 

Similar results are obtained with the simpler scalar model (jUj) - fl25|) . In this case, the 
continuum displacement vector of an edge dislocation at zero applied stress is 

U(x,y) = ^-9 (x - x , J J[)) ^ 



2tt V VA + 271 

where 6(1,0) is the Burgers vector and 9(x,y) = arctan(y/x). Choosing the singularity 
point (xq, yo) as explained above, we obtain configurations similar to those in Figures [2] and 
[3] except that there is no displacement in the vertical direction. To see the effect of a shear 
stress applied in the x direction, we select as initial and boundary condition 

B. Edge dislocation dipoles 

An edge dislocation dipole is formed by two edge dislocations with Burgers vectors in 
opposite directions. Depending on how we place the cores of these dislocations, different 
dipole configurations result. Let E(x, y) be the displacement vector (|26|) corresponding to 
the edge dislocation we have considered before. The static configuration corresponding to a 
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dipole is either a vacancy (Fig. @|, a divacancy (Fig. [5^) or a Stone- Wales defect (Fig. [5b). 
To obtain these different dipole cores, we use the following initial and boundary conditions 
at zero stress: 

• Vacancy E(x — x , y — y — 1/2) — E(x — x , y — y ). I = a/y/3 is the hexagon side in 
terms of the lattice constant a. 

• Divacancy E(a; -x ,y-y -l)- E(x -x ,y- y ). 

• Stone- Wales: E(x — x — a, y — y ) — E(x — x ,y — y ). 

We have set x = —0.25 a, y = —0.4/ to draw Figures H] and [51 

The vacancy represented in Fig. |4^b) is similar to that proposed in Ref. |24j using tight- 
binding calculations in graphite surfaces. The initial configuration represented in Fig. H^a) 
(asymmetric vacancy) evolves toward the dynamically stable symmetric vacancy represented 
in Fig. H](b) for overdamped dynamics. For conservative dynamics and zero initial velocity, 
the initial asymmetric vacancy evolves towards a stable oscillation about the symmetric 
vacancy. The symmetric vacancy has the threefold symmetry observed in experiments. In 
a recent paper, Telling et al.— propose that asymmetric vacancies are stable single vacancy 
defects in graphite sheets provided atoms are allowed to be displaced from the plane (see 



Figure 1 in Ref. 125). They further propose that, at room temperature, the displaced atom 
rotates around the vacancy center which would also explain the threefold symmetry observed 
in experiments (and possessed by the symmetric vacancy). For the initial conditions we have 
considered with either conservative or overdamped dynamics and at zero temperature, we 
have not found stable configurations resembling Fig. H](a) or oscillations between similarly 
asymmetric configurations with different orientations. Thus the asymmetric vacancy seems 
to be unstable for both overdamped and conservative dynamics in our model^. If this is 
indeed the case, this configuration is a saddle point with the conservative dynamics. Even 
if we allow bending modes of the graphene sheet, the instability of this saddle is likely to 
persist for small bending modulus. Future work will be devoted to study this problem and 
whether our stability results for different defects change if ripples in graphene are allowed. 

Both vacancies and divacancies are dynamically stable in an unstressed lattice. A suf- 
ficiently large applied stress along their glide direction (in both cases, the critical stress is 
about 0.12/i for a = 0.2) splits these dipoles, thereby originating two edge dislocations with 
opposite Burgers vectors that move in opposite directions. 
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Another possible core of an edge dislocation dipole is a 5-7-7-5 SW defect, as in Fig. [3(b). 
When introduced as initial condition of Eqs. ffT6l) - ffTTj) periodized with ffl9l) - (1211) . the SW 
configuration is likely to be unstable under zero applied stress. In fact, this configuration 
corresponds to two identical edge dislocations that have opposite Burgers vectors and share 
the same glide line. The dislocations comprising the dipole attract each other and are 
annihilated, leaving an undistorted lattice as the final configuration. Thus SW defects seem 
to be dynamically unstable^ except if we add external forces that produce the necessary 
bond rotation and stabilize them. If a shear stress is applied in their glide direction, these 
defects either continue destroying themselves or, for large enough applied stress (0.15/^ for 
a = 0.2), are split in their two component heptagon-pentagon defects that move in opposite 
directions as shown in Fig. [HI Note that stability of SW defects may be very different in 
small-radius nanotubes which, unlike graphene, do not have edges on their lateral surface 
and have a large curvature. 

If we shift the 5-7 defects comprising a SW in such a way that their glide lines are not 
the same, then we can obtain a configuration which is stable at zero applied stress. The 
simplest such case is a 7-5-5-7 defect in which the two pentagons share one side. We have 
checked that this defect is dynamically stable^ and that it splits in its two component edge 
dislocations when a shear stress exceeding 0.09 fi is applied to the lattice for a = 0.2. Note 
that in all cases (vacancies, divacancies, 5-7-7-5 and 7-5-5-7 defects) the critical shear stress 
needed to split the dipole is larger than the Peierls stress for an edge dislocation. The reason 
is that splitting a dipole requires overcoming the resistance of the lattice to motion (Peierls 
stress) and the attraction experienced by two edge dislocations of opposite Burgers vectors. 
The latter depends on the defect configuration which thus determines the critical dipole 
splitting shear stress. 

Instead of a dislocation dipole, our initial configuration may be a dislocation loop, in 
which two edge dislocations with opposite Burgers vectors are displaced vertically by one 
hexagon side: E(x — xq — a, y — y ) — E(x — xo,y — yo — I) (/ = a/^/3 is the length of 
the hexagon side). In principle, the dislocation loop could evolve to an inverse SW defect 
(7-5-5-7). Instead, this initial configuration seems to evolve towards a single octagon. If we 
displace the edge dislocations vertically by 1/2, E(x — x — a, y — yo) — E(x — x , y — y — l/2), 
the resulting dislocation loop seems to evolve towards a single heptagon defect. 
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C. Energetics 



We can have an idea about the energy associated to each defect by using the energy of 
the scalar model ffTSl) (measured with respect to the stress- free undistorted lattice), 

8 = (K«i) ~ M H] 2 + H n 2) - u(n)] 2 + [u{n 3 ) - u(n)f) 

n 

+ ~{\ + ») ([u(n 6 )-u(n)] 2 + [u(n 7 )-u(n)] 2 )}, (28) 

where we sum over all points in the hexagonal lattice, belonging to sublattices A or B in 
Fig. [7J, with neighbors given by d3J) or (fl2|) . respectively It can be seen that Eq. (TIB"]) is 
equivalent to pa 2 d 2 u/dt 2 = —d8/du. For a lattice with 16 x 16 lattice spacings, the energies 
associated to the defects we have described are 0.658 eV for the apparently unstable 5-7-7-5 
SW defect, 2.283 eV for the octagon with a dangling bond (shuffle dislocation), 4.917 eV for 
the 5-7 defect, and in the case of the dislocation dipoles: 8.825 eV for the vacancy, 12.074 
eV for the 5-8-5 divacancy and 9.483 eV for the stable 7-5-5-7 defectr^22. These values are 
similar to those found by Ewels et at^- for the activation barrier to form a glide dislocation 
dipole (8.99 eV) and a shuffle dislocation (2.29 eV) in graphene; cf. their Figure 3. Except for 
the SW defect, all other dislocation and dislocation dipole cores are stable and are obtained 
by dynamical evolution using the governing equations of the model from the class of initial 
conditions we mentioned above. 

VI. CONCLUSIONS 

We have studied edge dislocations and dislocation dipoles in planar graphene at zero 
temperature by means of periodized discrete elasticity models which seamlessly match the 
elastic field of dislocations and dipoles as the distance from their core increases. The cores 
of edge dislocations may be the well-known pentagon-heptagon defects of Fig. [2] (glide dis- 
locations) or octagons with a dangling bond (shuffle dislocations^) as in Fig. [31 depending 
on how we choose the initial configuration. Similarly, different cores are possible for edge 
dislocation dipoles: vacancies, 5-8-5 divacancies, Stone- Wales defects and 7-5-5-7 defects. Of 
these possible cores, symmetric vacancies, divacancies and 7-5-5-7 defects are dynamically 
stable whereas asymmetric vacancies and 5-7-7-5 SW defects are likely to be unstable. Our 
results show that regularizing linear elasticity near dislocation cores by periodized discrete 
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elasticity is a good alternative to computationally intensive atomistic simulations provided 
defects are sparse. 
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